//setup
   clear all
   set matsize 4000
   args SEED REPS TYPE
   capture cd "D:/DATA.IDB/Dropbox/1.Research/2.Methods/BDM/BDM2012/code/"
     if _rc==0 global dir = "D:/DATA.IDB/Dropbox/1.Research/2.Methods/BDM/BDM2012/code/"
   capture cd "D:/Dropbox/My Dropbox/1.Research/2.Methods/BDM/BDM2012/code/"
     if _rc==0 global dir = "D:/\Dropbox/My Dropbox/1.Research/2.Methods/BDM/BDM2012/code/"
   capture cd "/accounts/fac/jmccrary/research/BDM2012/code"
     if _rc==0 global dir = "/accounts/fac/jmccrary/research/BDM2012/code"
   cd "$dir"

//load library
   mata: mata set matastrict off
   run "ipw.ado"      //IPW1&IPW2
   run "gpe.ado"      //GPE
   run "llrmatch.ado" //LLR match
   run "mmatch.ado"   //NN&BCM based on X or pX
   run "frolichdgp.ado"  //Froelich DGP
   run "nswdgp.ado"      //NSW DGP
   run "msdgp.ado"      //DGP for misspecified estimators
   
//this program computes all estimators
//and it is called by all the simulation routines
capture program drop estimators
program define estimators, rclass
{ 
    syntax varlist(max=2), pscorex(string) matchx(string) biasadj(string) name(string)
    tokenize `varlist'
	local Y = "`1'"
	local T = "`2'"
    capture logit `T' `pscorex', iterate(50) asis
    if _rc==0{
	predict phat
  //IPW1
    ipw `Y' `T', pscore(phat) estimator(IPW1)
    return scalar ipw1_`name'= r(ipw1)
  //IPW2
    ipw `Y' `T', pscore(phat) estimator(IPW2)
    return scalar ipw2_`name'= r(ipw2)
  //GPE
    capture IPT_TOT `Y' `T' `pscorex', gen(That_IPT W_IPT That_logit W_logit)
    if _rc==0 return scalar gpe_`name'  = _b[TOT_`Y']
    if _rc!=0 return scalar gpe_`name'  = -9999999999
  //LLR match
    llrmatch `Y' `T', pscore(phat) kernel(epan)
    return scalar llr_`name'= r(llr)
  //NN & BCM match
    forvalues nn = 1(1)4{
      //match on phat
       mmatch `Y' `T' phat, m(`nn') biasadj(`biasadj') tc(att) kweight(KM)
       return scalar nn`nn'_pX_`name'= r(NN)
       return scalar bcm`nn'_pX_`name'= r(BCM)
      //match on X
       mmatch `Y' `T' `matchx', m(`nn') biasadj(`biasadj') tc(att) kweight(KM)
       return scalar nn`nn'_X_`name'= r(NN)
       return scalar bcm`nn'_X_`name'= r(BCM)
	 }  
   cap drop phat That_IPT W_IPT That_logit W_logit 
   }
   ereturn clear
}
end

//This program does Frolich sims (using a logit)
capture program drop simfrolich
program define simfrolich, rclass
{
   syntax, design(integer) type(string) [var(real 0.1)]
   frolichlogitdgp, n(100) design(`design') curves(1 2 3 4 5 6) var(`var')
   local pscorex = "Xk"
   // Match on X
   local matchx = "X"
   forvalues i =1(1)6{
      local biasadj= "m`i'"
   	  estimators y`i' T, pscorex(`pscorex') matchx(`matchx') biasadj(`biasadj') name(y`i')
	  return add
   }
   ereturn clear
}
end

//This program does Frolich sims (using a logit)
capture program drop simfrolich_est
program define simfrolich_est, rclass
{
   syntax, design(integer) type(string) [var(real 0.1)]
   frolichlogitdgp, n(100) design(`design') curves(1 2 3 4 5 6) var(`var')
   if `design'!=3  local pscorex = "C1 C2 C3 S1 S2 S3"
   if `design'==3 local pscorex = "C1 C2 C3 C4 C5 S1 S2 S3 S4 S5"
   // Match on X
   local matchx = "X"
   forvalues i =1(1)6{
      local biasadj = "${reg`i'}"
   	  estimators y`i' T, pscorex(`pscorex') matchx(`matchx') biasadj(`biasadj') name(y`i')
	  return add
   }
   ereturn clear
}
end


//This program does NSW sims
capture program drop simnsw
program define simnsw, rclass
{
   syntax, k(integer) type(string)
   nswdgp, k(`k') replicas(1)
   local matchx  = "age educ earn74 earn75 married unempl74 unempl75 dropout"
   local pscorex = "age educ earn74 earn75 earn74_earn75 married unempl74 unempl75 dropout unempl74_unempl75 earn74sq earn75sq"
   local biasadj = "age educ earn74 earn75 earn74_earn75 married unempl74 unempl75 dropout unempl74_unempl75 earn74sq earn75sq"
   estimators y1 T, pscorex(`pscorex') matchx(`matchx') biasadj(`biasadj') name(y1)
   return add
   ereturn clear
}
end

//This program does misspec sims
//It first gets the data (curve 4, design 1, fourier flavor), then computes all estimartors
//looping over different specifications of the pscore and bias-adj model
capture program drop simmisspec
program define simmisspec, rclass
{
  syntax, wiggle(real)
  msdgp, n(400) wiggle(`wiggle')
  forvalues i =1(1)4{
  local lin   = "X1 X2 X3 X4"
  local inter = "X1X2 X1X3 X1X4 X2X3 X2X4 X3X4"
    if `i'==1 local pscorex = "g" 
    if `i'==2 local pscorex = "`inter'"
    if `i'==3 local pscorex = "`lin'"
    if `i'==4 local pscorex = "`lin' `inter'"
  forvalues j =1(1)4{
    if `j'==1 local biasadj = "m" 
    if `j'==2 local biasadj = "`inter'"
    if `j'==3 local biasadj = "`lin'"
    if `j'==4 local biasadj = "`lin' `inter'"
    estimators y1 T, pscorex(`pscorex') matchx(`lin') biasadj(`biasadj') name(y1_`i'_`j')
    return add
   }
   }
   ereturn clear
}
end

//Simulations for tables 1, A1a, A1b: Frolich (Simple and well specified models)
if "`TYPE'"=="frolich" | "`TYPE'"=="frolich001"{
    local start = "$S_TIME"
    set seed 1234567`SEED'
	cap log close
	log using "$dir/sim`TYPE'_`SEED'.log", replace text
	if "`TYPE'"=="frolich"     local var = 0.1
	if "`TYPE'"=="frolich001"  local var = 0.01
	forvalues d = 1(1)5{
	   simulate, reps(`REPS'): simfrolich, design(`d') type(`TYPE') var(`var')
	   gen design = `d'
	   if `d'==1 save "$dir/sim`TYPE'_`SEED'.dta", replace
	   if `d'> 1{
				 append using "$dir/sim`TYPE'_`SEED'.dta"
				 save "$dir/sim`TYPE'_`SEED'.dta", replace
	   }
	}
	local end = "$S_TIME"
	di "Started at `start'. Finished at `end'."
}

//Simulations for tables 1, A1c, A1d: Frolich (Simple and well specified models)
if "`TYPE'"=="frolich_est" | "`TYPE'"=="frolich_est001"{
    local start = "$S_TIME"
    set seed 1234567`SEED'
	cap log close
	log using "$dir/sim`TYPE'_`SEED'.log", replace text
	if "`TYPE'"=="frolich_est"     local var = 0.1
	if "`TYPE'"=="frolich_est001"  local var = 0.01
	forvalues d = 1(1)5{
	   simulate, reps(`REPS'): simfrolich_est, design(`d') type(`TYPE') var(`var')
	   gen design = `d'
	   if `d'==1 save "$dir/sim`TYPE'_`SEED'.dta", replace
	   if `d'> 1{
				 append using "$dir/sim`TYPE'_`SEED'.dta"
				 save "$dir/sim`TYPE'_`SEED'.dta", replace
	   }
	}
	local end = "$S_TIME"
	di "Started at `start'. Finished at `end'."
}

	
//Simulations for tables 2, A2: NSW (Well specified models with different overlap parameters)
if "`TYPE'"=="nsw"{
	local start = "$S_TIME"
	set seed 12345678`SEED'
	cap log close
	log using "$dir/sim`TYPE'_`SEED'.log", replace text
	//Bad Overlap
	simulate, reps(`REPS'): simnsw, k(1) type(`TYPE')
	gen overlap = 1
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	//Good Overlap
	simulate, reps(`REPS'): simnsw, k(5) type(`TYPE')
	gen overlap = 5
	append using "$dir/sim`TYPE'_`SEED'.dta"
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	local end = "$S_TIME"
	di "Started at `start'. Finished at `end'."
 }

 
//Simulations for tables 4, A3: Misspecified models
if "`TYPE'"=="misspec"{
    local start = "$S_TIME"
	set seed 1234567`SEED'
	cap log close
	log using "$dir/sim`TYPE'_`SEED'.log", replace text
	//Linear DGP
    simulate, reps(`REPS'): simmisspec, wiggle(0)
	gen wiggle = 0
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	//Interactions DGP
    simulate, reps(`REPS'): simmisspec, wiggle(1)
	gen wiggle = 1
	append using "$dir/sim`TYPE'_`SEED'.dta"
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	//Linear+Interactions DGP
    simulate, reps(`REPS'): simmisspec, wiggle(0.5)
	gen wiggle = 0.5
	append using "$dir/sim`TYPE'_`SEED'.dta"
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	//Linear selection eq +Interactions ouctome eq DGP
    simulate, reps(`REPS'): simmisspec, wiggle(10)
	gen wiggle = 10
	save "$dir/sim`TYPE'_`SEED'.dta", replace
	//Interaction selection eq + linear ouctome eq DGP
    simulate, reps(`REPS'): simmisspec, wiggle(20)
	gen wiggle = 20
	append using "$dir/sim`TYPE'_`SEED'.dta"
	save "$dir/sim`TYPE'_`SEED'.dta", replace
}


/// ------> To run the simulations open SIMULATIONS <--------

